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We present MadDM v.2.0, a numerical tool for dark matter physics in a generic model. 
This version is the next step towards the development of a fully automated framework for 
dark matter searches at the interface of collider physics, astro-physics and cosmology. It 
extends the capabilities of v.1.0 to perform calculations relevant to the direct detection of 
dark matter. These include calculations of spin-independent/spin-dependent nucleon scat¬ 
tering cross sections and nuclear recoil rates (as a function of both energy and angle), as 
well as a simplified functionality to compare the model points with existing constraints. The 
functionality of MadDM v.2.0 incorporates a large selection of dark matter detector mate¬ 
rials and sizes, and simulates detector effects on the nuclear recoil signals. We validate the 
code in a wide range of dark matter models by comparing results from MadDM v.2.0 to the 
existing tools and literature. 
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I. INTRODUCTION 

While there are many intriguing, indirect hints of the existence of dark matter (DM), the evi¬ 
dence of irrefutable signals and the true nature of DM remain elusive (see Ref. [I] and references 
therein for recent progress in DM searches). Various direct detection experiments are currently on¬ 
going, and many are under development. The results from direction detection experiments have so 
far produced stringent constraints on the DM scattering cross section off atomic nuclei, especially 
for the DM in the mass range of the electroweak (EW) scale. More recently, the LHC analyses 
of jet(s) -|- missing energy and WjZjt -|- missing energy channels, and their variants, provided 
complementary information to underground fixed target experiments. The complementarity was 
particularly evident in the improved sensitivity for a relatively light DM in simplified model sce¬ 
narios [2]. Despite a significant effort, none of the individual experiments have provided a strong 
indication on the mass scale of dark matter particles. It is therefore essential to search for evidence 
of dark matter in many different ways: by directly producing DM candidate at particle colliders, 
as well as by detecting it in our galaxy and beyond [3]. 

Searching for DM at the interface of collider physics, astroparticle physics and cosmology is 
perhaps the most promising path to discovery and discrimination of different DM models, which 
requires substantial amount of computational power and relevant tools. While many collider event 
generators are publicly available |1], the number of numerical tools which are able to calculate 
signals of dark matter at the interface of collider physics and astroparticle physics remains very 
limited urn- MadDM [7] emerged as an attempt to build such a tool. The goal of the MadDM 
project is to form a simple, user friendly and all-in-one DM phenomenology framework. Our 
ambition is to allow both experimentalists and theorists to calculate accurate signatures of generic 
DM models at colliders, in our galaxy and in the early universe with only a few clicks of a button. 
We also aim to provide flexibility such that users easily replace the existing modules with their 
own implementation or link MadDM functionality with other existing tools. 

The user friendly architecture of the MadGraph5_aMC@NLO UM package provides an ideal 
framework for the development of MadDM. Numerous MadGraph additions and extensions offer 
many interesting directions that development of MadDM can go into. Version 1.0 of MadDM 
focused on computation of DM relic density in a generic model [7]. Here we discuss the next 
step in the development of MadDM: direct detection of galactic DM. In addition to the ability to 
calculate DM-nucleon cross sections, MadDM v.2.0 goes a step further and provides the directional 
information of nuclear recoil, as well as the ability to simulate the effects of detector systematics 
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on recoil distributions. To the best of our knowledge, MadDM is the first publicly available code 
with such a capability. The code is able to calculate double differential recoil rates, as a function 
of recoil angle and energy, as well as the energy distribution and angular distribution of nuclear 
recoil. 


Directional information of DM scattering can be of great importance for the next generation 
dark matter direct detection experiments. In case a DM signal is observed, directional information 
about nuclear recoil rates could lead to more precise measurements of intrinsic particle properties 
of DM |lUj as well as astrophysical information such as the DM velocity distribution jllL 112] . 
Conversely, if next generation direct detection experiments lead only to more stringent limits on 
DM mass and scattering cross section, directional detection could become crucial in overcoming 
the irreducible neutrino background m- 


Finally, MadDM v.2.0 incorporates a simplified functionality for testing model points against 
experimental constraints, whereby the user can choose the range of DM relic densities and upper 
bounds on DM-nucleon scattering cross sections for which the model is consistent with the exper¬ 
iments. The code will automatically compare the output of the calculation to the user specified 
constraints and determine whether the model point is allowed or not. 


New version of MadDM follows the already established structure of MadDM v.1.0. A Python 
module generates relevant amplitudes for relic density and direct detection (with ALOHA [14] ). 
while a FORTRAN module handles the heavy numerical calculations. As before, MadDM is compat¬ 
ible with UFO (Universal FeynRules Output [15]) model which contains a dark matter candidate 
and can be easily linked to any width or mass spectrum generator which can produce a Les Houches 
formatted parameter card miiiTi. 


The following sections describe the new functionality of the MadDM code for the direct detection 
of DM, while we refer the reader to Ref. [7| for more information on calculation of DM relic 
density. Section |TI| describes the theoretical background of dark matter direct detection, including 


directionality. A description of the new routines in the MadDM code is presented in Section III 


Finally, we show several validations and example calculations using the MadDM v.2.0 code in 


Section IV We reserve the appendix for a brief description of the MadDM code structure. 
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II. REVIEW OF DIRECT DETECTION AND ITS IMPLEMENTATION IN MADDM 

A. DM-nucleon Elastic Cross Section 

The possibility that DM could weakly interact with ordinary matter has led to a development 
of a myriad of experiments which search for signatures of galactic DM scattering off nuclei |18H22] . 
Being sensitive to the nuclear recoil of the scattering DM-nucleon event in underground detectors, 
experiments can constrain model parameters involved in the DM-nucleon (nucleus) elastic scatter¬ 
ing cross section, while in the case of a discovery, the same information may be used to determine 
properties of DM particles. The abundance of independent experiments has so far led to impor¬ 
tant limits on the DM-nucleus cross-section, among which the most recent (and most stringent) 
one comes from the LUX collaboration [I8l[23]. In this section, we give an overview of the theo¬ 
retical formalism of DM-nucleus elastic scattering in order to allow the reader to understand the 
procedures MadDM performs in computing DM direct detection rates [24H26| . 

The MadDM v.2.0 code incorporates the following definition of the DM-nucleus 
spin-independent (SI) scattering cross section: 

= —A ' fp + ' /n]^ , (1) 

TT 

where = m^^ruA DM-nucleus reduced mass and fp and fn are proton/neutron spin- 

independent form factors respectively [IIl[l3l[27j. We use Z to denote the number of protons in 
a nucleus and A to denote the number of protons and neutrons. Similarly, for spin-dependent 
(SD) interactions, we use the definition: 

'^SD ^ \Jp T Jn) ) (2) 

where fp and are the proton/neutron spin-dependent form factors respectively, and Ja is the 
spin of the nucleus. Particle physics enters the calculation of ctsd/si the form factors /at and 
f'j^ where N = p,n. 

The spin-independent quantity for a scalar current, /jy is defined as 



where aq are the coefficients of the low energy matrix elements defined by 


{M) ixq Xq) = Oiq {ipgll^q) ■ 


( 4 ) 
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The quantities are nucleon form factors related to the low energy elastic scattering matrix 

elements as follows 

{N I'lljg'ijjql N) = — {ior u,d,s), ( 5 ) 

Ulq 

{N \'>pqiljq\ N) = (for c, b, t). (6) 

The values for nucleon form factors are typically extracted from data or from chiral perturbation 
theory [28] and depend on the type of current which defines the DM-nucleon interaction. The 
default values for scalar spin-independent form factors of light quarks used in MadDM are: 

= 0.0153, = 0.0110, 

= 0.0191, = 0.0273, (7) 

/f = 0.0447, = 0.0447, 

while the gluon scalar form factor is 

fN' = '^-'EfN- ( 8 ) 

Q 

The vector current interaction is characterized by a simpler set of form factors. Since the 
vector current is conserved, the DM-nucleon form factor can be obtained from the sum of the 
currents of the valence quarks, i.e.: 

fN= Y. ( 9 ) 

q=u,d 

where = /^'^ = 2 and = 1. 


For spin-dependent interactions, the following equation defines the proton/neutron form factors 
I'n in Eq. 


/at— Yj 

q=u^d^s 

where the quantities are nucleon form factors, applies for both axial-vector and tensor 
currents. In MadDM, we use the following default values for the axial-vector interaction: 


0.842, 

-0.427, (11) 

A^^" = A;^^" = -0.085, 
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while we use another set of coefficients for the tensor interaction: 

= A™ = 0.84, 

A™ = A^“=-0.23, (12) 

AJ" = = -0.046. 

The nucleon form factors for the axial-vector current, Ap^“ and have been provided by 

HERMES and COMPASS experiments |29l [30] while tensor coefficients come from lattice calcula¬ 
tions [3T] . 

The coefficients aq can be extracted from the full matrix elements for XQ scattering via Fiertz 
transformations. While this is relatively straightforward to do analytically, numerical implementa¬ 
tions of Fiertz transformations are non-trivial. In the following section, we discuss the numerical 
method MadDM uses to compute the low energy coefficients aq from the full matrix elements in 
the XQ scattering. 


B. Projection Operator Method for Extraction of Low Energy Coefficients 


The procedure for computing the low energy coefficients for DM-nucleon scattering in MadDM 
is similar to the procedure implemented in MicrOMEGAs [2^ At the quark level, the input 
quark-DM interacting Lagrangian at low energy can be re-written as a set of effective operators as 
follows: 

= (13) 

i 

where the operators Oq standing for xQ scattering are defined in Table [I| and the index i runs over 
all the relevant contributions. The list of operators can be separated in odd and even operators 
under the interchange of quarks and anti-quarks as follows: 


-^^q2=o — (-^s/ + ^'si) + (-^Id + ^°Sd) 


q,s 


SIraSI 


^SD^nSD 




q,s 


(14) 


where q = u,d, s, c, b, t and s stands for even (e) or odd (o) operator. The Fiertz transformation 
which gives aq^s coefficients projects the matrix element for xQ XQ onto the set of effective 
operators. Numerically, this is equivalent to taking the interference term of the matrix element 


^ For improved treatments on QCD effect in dark matter scattering off nucleon, see Ref. |32| and Ref. [33] for recent 
development. 
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DM spin 

Even 

Odd 



scalar current 

vector current 

SI 

0 







1/2 

V'xV'xV'gV’g 



1 

2M^A^^A'^'tpq'll)q 

i {A*^d^A^a - A“5^A*„) 



axial-vector current 

tensor current 

SD 

1/2 



1 

x/6 (^daA*^pA^^ - A*^pdaA^^'^ 



Table I. List of effective operators taken from Ref. [26], implemented into MadDM. S, ip-^ and correspond 
to scalar, fermion and vector DM fields, respectively. 


for XQ XQ scattering with the even effective operator matrix element for the same process. The 


projection over the effective operators will select either SI or SD contributions since Eq. (13) is 
written in an orthogonal basis. 

The even and odd coefficients can be further separated by considering scattering off of both 
quarks and anti-quarks. The sum and difference of the even and odd coefficients can be written as: 


a, 


q,e 


T OLq,o — 


\M^* ■ 




(15) 

(16) 


where the lAlg^p in the denominator takes into account the fact that the effective operators are 


not properly normalized. In Eq. (16) we used the property that the odd coefficient changes sign 
under the interchange of g —?■ g. Note that the procedure of projection operators for extracting aq^s 
is completely generic although the list of effective operators is non-exhaustive. Indeed, operators 
such as ^ 75 V’ are neglected since they are suppressed in the zero momentum limit. 


C. Differential Nucleus Recoil Rate 

The direct detection module of MadDM allows the user to calculate quantities relevant for 
DM-nucleus scattering beyond the spin-independent and spin-dependent cross section. One is, 
for instance, often interested in studies of nucleus recoil rates, where traditionally the quantity of 
interest has been dR/dE, the differential nuclear recoil rate as a function of recoil energy. However, 
valuable information is contained in the angular distribution of the recoil rate dR/dVt. And despite 
the fact that no currently operating experiment has the ability to efficiently measure directionality 
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of the nuclear recoil with high statistics, there are few ongoing efforts in this direction as well 
as concrete proposals for future experiments [ini [la EiHST]. Hence, we deemed it important to 
include the ability to calculate the angular rates into the MadDM code. 

One reason why angular information could be important is that low energy recoil events recorded 
by direct detection experiments lie well within the range dominated by cosmogenic backgrounds 
and those resulting from radioactivity. A strong signature of dark matter may come from direc¬ 
tional detection, which exhibits a very large forward-backward asymmetry [Afb) in the angular 
distribution of nuclear recoil events. A large Aps would be difficult to mimic with any known 
backgrounds as one expects the angular distribution of such processes to be roughly uniform 

The fact that background processes are (typically) spatially isotropic could be used to overcome 
another impeding problem. As the direct detection bounds on the DM-nucleon cross-sections get 
lower, neutrino backgrounds will present an obstacle [38] and directional detection could provide 
a way to circumvent the so called “neutrino floor” |39) . 

In case a signal is ever observed, the study of directionality in DM-nucleus recoils will be of 
foremost interest. The angular distribution of the recoil rates could provide useful information 
about the astrophysical aspects of the DM galactic halo prohle nil SO], as well as information on 
the anisotropy in the DM halo velocity distribution and prospects for more accurate measurements 
of the DM mass and interaction cross-section m- 

In this section we present a detailed overview of the theoretical background of differential recoil 
rates for DM detection. We begin by considering our solar system moving through the Galactic 
“DM-wind” in the direction of the Cygnus X — 2 constellation as illustrated in Fig. A DM 
particle of mass m^, incident at velocity v = u(sin a cos /? x -|- sin a sin /3 y -|- cos a z) in the detector, 
elastically scatters off a target nucleus of mass mjv, stationary in the detector, and causes it to 
recoil with velocity u = M(sin0cos(/)x -|- sin0sin(/>y -|- cos0z) and momentum q in the direction 
{6,4>). The event rate per unit recoil energy and unit recoil solid angle is 




2/30 


da 


-vf{v)d^v, 


(17) 


dERdn(^g^^) J dq^dn^^g^^) 

where pQ is the dark matter density in our local galaxy and /(v) is the velocity distribution of the 
DM in the galactic halo. The differential cross-section 


da 


Sirp^v 


F‘^{q)6[v cos 9 



(18) 


is obtained from the two body scattering phase space distribution [T2|. Here, p = m^rriN/{my. + 
mjxf) is the reduced mass of the DM-nucleus system, q = y/2 Er is the recoil momentum 
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Figure 1. Geometry of DM scattering off a target nucleus in the detector. The DM is incident at an angle 
(a, /3) relative to the z axis. The nucleus recoils in the {0, (p) direction. Note that (3 and p angles lie in 
the X, y plane. The thin and dashed thin lines serve to illustrate the projection of the incoming and recoil 
vectors onto the z axis and the xy plane. 


with E[i the recoil energy, angle 0 is the recoil angle which determines the direction between the 
recoiling nucleus and the initial DM trajectory, while CTy-jy is the DM-nucleus scattering cross- 
section, weighted by the nuclear elastic scattering form factor F{q). Assuming the nucleus is a 
sphere with uniform density, the form-factor is the Fourier transform of the nuclear density and 


includes all the relevant nuclear effects. For a detailed discussion on F{q) see Section IID The 
DM-nucleus cross-section is summed over contributions from the spin-independent {(^si) and spin- 
dependent {(tsd) cross-sections respectively, as defined in Eq. Q and Eq. Q respectively. Eq. 0 
becomes 


d?R 


dFdQf 


TrArmi J 2^ 


(19) 


‘( 0 , 0 ) ' "-x 

where r = AmNiTix /{tun + Nq = 6.022 x 10^® kg“^ is Avogadro’s number, = 0.932 A 

is the target mass and A the atomic mass number (the factor of 0.932 is the value of atomic mass 
units (AMU) in GeV). This form of the double differential recoil rate is particularly useful as the 
quantity 


f{VqA)= / 5iy^-Vq)fiy)d^V, 


( 20 ) 
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is easily identified as the three dimensional Radon transforn:^ where q is the recoil momentum 
direction, v is the velocity of the DM particle and Vq = -^ is the minimum velocity required for the 
DM to impart a recoil momentum q on the nucleus. f{vq, q) represents the velocity distribution for 
a detector stationary in the galactic frame. We consider an observer moving through the Galaxy 
with velocity Viab in the galactic frame. Viab is related to the DM velocity in the lab frame (viab) 
and the DM velocity in the galactic frame (vgai) by the Galilean transformation viab = Vgai — Vjab- 
The transformation takes into account the motion of a detector relative to the DM particles as 
seen in the lab frame. By the properties of the Radon transform for a pure translation Viab, we 
obtain 


fuh{Vq, q) = fgalivq + Viab-A, 4)- 


( 21 ) 


We assume a Maxwell-Boltzmann velocity distribution which is truncated at the DM escape velocity 
Vesc such that 

1 


fuBiv) = 






( 22 ) 


for V < Vesc and fuBiv) = 0 otherwise. In Eq. (22), vq is the most probable speed of the DM in 
the halo with typical values of 220 ~ 250 km/s and kesc is a normalization factor which is obtained 
by integrating the velocity distribution in the galactic frame from 0 to Vesc- For Viab = where 


ve is the velocity of the Earth with respect to the galactic frame, Eq. (21) becomes 

{Vq + VE.q)^ 


fuBiVq+^E-qA) = 


1 


kesc^/^VQ 


exp 


exp 


(23) 


Gombining Eq. (19) and Eq. (23) we obtain an expression for a double differential recoil rate in 

{VE COsO - Vrmnf' 


the laboratory frame: 

(PR _ 2Nopoa^N F\Er) 


dEjid cos 0 Arvo ke 


exp 


exp 


(24) 


^0 / \ 

where ve • q = —vecosO and Vq = Umin = vq\/Er/Eqv, with Eq = l/2m^UQ being the most 
probable kinetic energy of the DM. The velocity of the Earth ue, with respect to the galactic 
frame is calculated in the appendix of Ref. [JT] and includes the motion of the Earth with respect 
to the Sun, the proper motion of the Sun and the motion of the solar system with respect to the 


galactic center. Eq. (24) assumes azimuthal symmetry of the DM velocity profile as illustrated in 
Fig.[Tj 


For more information about the Radon transforms as applied to directional detection, we refer the reader to 

Ref. [12]. 
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D. Nuclear Form Factors and Detector Resolution 


In order to compute the recoil rates, we take into acconnt both spin-dependent (axial vector) and 
spin-independent (scalar) DM-nucleon interactions. The contributions of different types of inter¬ 
actions show up in the DM-nucleus cross-section and in the form-factor. For the spin-independent 
case, we consider the Helm form factor as described in Ref. mi: 

FsM = , 25 ) 

with q = 6.92 x fm“^ the recoil momentum, + 7a? — fm the 


effective radius of the nucleus, c = (1.23A^/^ — 0.60) fm, a = 0.52 fm and s = 0.9 fm. Eq. (25) is 


in fact the Wood-Saxon form factor, with = (? — 5s^ and c = 1.2A^/^ fm. Furthermore for the 
spin-dependent case, we consider 

S{q) 


FiBiq) = 


5(0)’ 


(26) 


OO — + ^d.) + ’ 

ai = ^(A,2 + A,) +Aa,, 


where the recoil momentum dependent structure function is given by S{q) = aQSoo{q)+aoaiSoi{q) + 
a1Sii{q) and 5(0) is the structure function at zero momentum transfer. For the coefficients 

(27) 

(28) 

we use the nucleon form factors, Aj {i = u,d,s), as defined in Eq. 0- The nuclear spin 
structure functions Sij for different nuclei are approximated using htting functions. Eor Xenon, 
Iodine and Sodium we use the functions from the Nijmegen II fitting as calculated in Ref. |42] . 
Eor Germanium we use those calculated in Ref. [l3], while for Fluorine and Silicon we use the 
fitting functions found in Ref. [H] . All the information we use for both spin-independent and spin- 
dependent interactions are summarized in Table |ll] and Table m respectively. For the nuclear 
structure functions that are cnrrently not very well calculated (such as Carbon, Sulphur, Neon and 


^2 p2 


as described in Ref. 


Argon), we use the Gaussian distribntion SM = Si, (0)exp 
with Ra = 1.7AV3 - 0.28 - 0.78(AV3 _ 3 g ^(AV3 _ 3 . 8)2 -h 0.2) A < 100 and Ra = 1.5AV3 
for A > 100. We note that the Gaussian approximation breaks down for large recoil energies 


Finally we use Eq. (24) to simulate the scattering of dark matter with the nuclens in a detector. 


We take into account the energy and angular resolutions of a generic detector as applied to the 


® We were also not able to find the values in asterisk. For Neon, Argon and Snlphur, we set the orbital spin values 
to I to avoid any infinities in the calculation of the spin-dependent cross-sections as shown in Eq. ([^, while we 
set the default values of the magnetic moments for Argon and Sulphur to 0. 
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Target Material 


List of Stable Isotopes & Abundances 

Xenon 

A 

123.9, 125.9, 127.9, 128.9, 129.9, 130.9, 131.9, 133.9, 135.9 

Z = 54 

Abundance (%) 

0.09, 0.09, 1.92, 25.44, 4.08, 21.18, 26.89, 10.44, 8.87 

Germanium 

A 

69.92, 71.92, 72.92, 73.92, 75.92 

Z = 32 

Abundance (%) 

20.84, 27.54, 7.73, 36.28, 7.61 

Silicon 

A 

27.98, 28.98, 29.97 

Z = 14 

Abundance (%) 

92.23, 4.68, 3.087 

Neon 

A 

19.99, 20.99, 21.99 

Z = 10 

Abundance (%) 

90.48, 0.27, 9.25 

Argon 

A 

37.962, 39.962 

Z = 18 

Abundance (%) 

0.0632, 99.6 

Sodium 

A 

22.9 

Z = 11 

Abundance (%) 

100 

Iodine 

A 

126.9 

Z = 53 

Abundance (%) 

100 

Carbon 

A 

12.0, 13.0 

Z = 6 

Abundance (%) 

98.89, 1.11 

Fluorine 

A 

18.998 

Z = 9 

Abundance (%) 

100 

Sulphur 

A 

31.9, 32.97, 33.96, 35.96 

Z = 16 

Abundance (%) 

94.9, 0.76, 4.29, 0.02 


Table II. List of stable isotopes and their abundances for the different target materials used in the calculation 
of the spin-independent quantities. 

distribution of the scattering events. We assume a Gaussian resolution function to smear the 
distribution of recoil events in both energy and angle with standard deviations aE = AV^ (the 
energy resolution) and (constant angular resolution) as in Ref. [TO]. Default resolutions are 
given by A = 1 and ag = 30° with an option for users to implement their own detector resolutions. 

III. HOW TO USE MADDM V.2.0 

A. Running MadDM 

Addition of the direct detection code to the existing MadDM package follows the syntax and 
philosophy of the previous MadDM version. The code is split up into the Python and the FORTRAN 
modules, whereby diagram generation and the structure of folders and files is handled by Python, 
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Target Material 


{Sn) 


J 

Xenon 

Ref. [H] 

-0.272 

-0.009 

3 

2 

Germanium 

Ref. [43] 

0.378 

0.030 

9 

2 

Silicon 

Ref. HD 

0.13 

-0.002 

1 

2 

Neon 

Refs. pFl Hi] 

0.294 

0.02 

1 * 

2 

Argon 

Ref. [26] 

0* 

0* 

1 * 

2 

Sodium 

Ref. [H] 

0.0199 

0.2477 

3 

2 

Iodine 

Ref. [H] 

0.075 

0.309 

5 

2 

Carbon 

Ref. [26] 

-0.172 

-0.009 

1 

2 

Fluorine 

Refs, nuns] 

-0.0087 

0.4751 

1 

2 

Sulphur 

Ref. [26] 

0* 

0* 

1 * 

2 


Table III. List of structure functions, moments and angular momenta for different materials. The fitting 
polynomials for the structure functions can be found in the references cited. We were not able to find the 
magnetic moments for Argon and Sulphur thus we set them to be zero, thus MadDM does not currently 
compute the SD rates for them. We were also not able to find the orbital spin values for Argon, Neon and 
Sulphur, we set these to ^ in the code at the moment. 

and the numerical calculations are performed by the FORTRAN module. As in the previous ver¬ 
sion, the code does not require any pre-compilation. Simply placing the code within a Mad- 
Graph5_aMC@NLO folder suffices. Note however, that Numpy is now a pre-requisite for running 
the MadDMv.2.0 code. In order to use the plotting functionality, Matplotlib is also required. 

In the following we briefly describe how to use the MadDM code relevant for direct detection 
of dark matter, while we give detailed descriptions of the main functions/features of the MadDM 
v.2.0 package in the Appendix]^ 

As in version I.O, the main MadDM program can be executed via the maddm.py script. The 
interface will prompt the user to enter the name of the UFO dark matter model (TAB auto-complete 
included), a name of the projects, as well as offer a choice of available computations (relic density, 
DM-nucleon cross sections and recoil rates). The Python module of MadDM will then proceed to 
generate all the relevant scattering/annihilation diagrams and set up the directory structure of the 
user defined project in Projects/<projectname> folder. 

The user can choose to enter the DM candidate manually, or to allow MadDM to determine the 
dark matter candidate automatically (the code assumes that the lightest particle with the PDG 
code greater than 25 and zero width is the DM candidate). The code also offers the user to edit the 


For the reference on the relic density calculation see Ref. [7]. 
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default parameter card, and pre-set the model parameters before it searches for the DM candidate. 

Upon determining the DM candidate, MadDM will attempt to find any co-annihilation candi¬ 
dates. The user is prompted to enter the mass splitting within which the co-annihilation partners 
are deemed to contribute significantly to the calcnlation of relic density (the default is 10 %). 

Finally, the code will prompt the user to edit the niaddm_card. inc file, which contains all 
relevant parameters for the numerical calculations. These include, but are not limited to precision 
parameters for the relic density calcnlation, nucleon form factors for direct detection, choices of 
target materials and detector parameters. Upon execntion and output of the results, the code 
will offer the user the option to perform a parameter scan. Alternative to the nser interface in 
maddm.py the MadDM code can also be run within any user defined Python script. We provide 
such an example in the example.py file. 

B. UFO Model Conventions 

Now we proceed to discuss some important conventions on UFO files reqnired by the Mad- 
DMv.2.0 code. We implemented the projection operator formalism in MadDM v.2.0 via a recent 
feature of MadGraph5_aMC@NLO which allows the user to merge two UFO model files. We pro¬ 
vide UFO models with vertices of effective operators as a part of the MadDM code, so that an 
interference term between the matrix elements of the user defined model, and the effective oper¬ 
ators can be extracted. The vertices are then automatically added to the DM model provided by 
the user during the execution of the code. The merging procedure of MadGraph5_aMC@NLO is 
sensitive to the UFO convention and the user should take special care to make sure that the DM 
model used in MadDM is compatible with the effective operator models embedded in MadDM 
v . 2.(|3 We have generated our set of effective operators with FeynRules 2.0. We advise the user to 
test the merging procednre in MadGraph5_aMG@NLO before rnnning MadDM in order to avoid 
UFO conflicts. Using this UFO convention, the DM particle should be defined as 

DM = Particle(pdg_code = pdg code, 

spin = 1, 

mass = Param. DMMASS, 
width = Param. WDM, 

® The statement is relevant only for direct/directional detection. Relic density will proceed correctly without regard 
for the UFO version. 



16 


...) 

in case of a real scalar DM (for example). In this convention, it is compulsory to assign a 
PDG code to the DM particle as well as a mass and a width. Next, the mass and width 
block parameters have to match the following template: 

DMMASS = Parameter(name = ’DMMASS’, 

lhablock = ’MASS', 
lhacode = [ pdg code ]) 

WDM = Parameter(name = ’WDM’, 

lhablock = ’DECAY’, 
lhacode = [ pdg code ]) 

where the lhablock and lhacode parameters should be exactly like the template when replacing 
the PDG code by the PDG code of the DM particle. Note also that we cannot use mdm for the 
DM mass name due to conflict with internal MadDM parameters. For the same reason, if the 
user wants to define new coupling orders, the following ones should not be used ; SIEFFS, SIEFFF, 
SIEFFV, SDEFFF, SDEFFV. 

Finally, the user should make sure that direct detection diagrams have a QED order 
lower or equal to 2 since MadDM generates diagrams with the QED=2 flag. 

IV. VALIDATIONS 

We performed validations of the MadDM code in a diverse class of benchmark models. In each 
case, we verified that the code reproduces accurately the results from existing literature or public 
codes such as MicrOMEGAs. 

A. Validations of the DM-nucleon Cross Section Calculation 

1. Simplified DM Models 

The simplest benchmark model we consider for validation is a scalar DM which communicates to 
the SM via the Higgs boson. The model contains only one diagram contributing to direct detection 
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i.e. a t-channel Higgs exchange. The Higgs-DM coupling is included in the ^H^HS operator where 
S is the DM scalar field. It implies that the model leads only to spin-independent cross section. 

For the purpose of validating calculations of the spin-dependent cross sections, we consider a 
similar simplified model, with an axial-vector mediator which couples through fermionic DM and 
SM via the following Lagrangian: 

Xl'"l5XVti + 9q , (29) 

where x and are the DM and mediator fields, respectively. Fig. shows our results for 
proton/neutron DM cross section for the two simplified models against results obtained with Mi- 
crOMEGAs. In both cases, we find an excellent agreement between MadDM and MicrOMEGAs 
with the differences at < 1% level. We have also checked that the results from MadDM are con¬ 
sistent results from MicrOMEGAs in case of simplified models with different mediators (scalar, 
fermion and vector) and different DM candidate (scalar, fermion and vector), which are shown in 
Table m 


2. Minimal Universal Extra Dimensions 


After validating MadDM for a class of simplified models, we turn to more complex models. 
For this purpose, we chose to examine the results of MadDM for the Minimal Universal Extra 
Dimension model (MUED) against the results in existing literature. 

MUED is the simplest model containing a Kaluza-Klein (KK) dark matter candidate among 
extra dimensions theories. The lightest KK particle, the KK-photon ( 71 ) appears as a dark matter 
candidate in vanilla UED scenarios. At tree level, the inverse radius R~^ of the extra-dimension 
corresponds roughly to the mass of the massive fields at level one. The mass splitting between 
KK-quark (gi) and KK-photon: 




m. 


(30) 


'7i 


plays an important role as a free parameter in direct detection. The details and phenomenology 
of MUED will not be discussed here since it has been studied extensively in the past |46II49j and 
we will proceed directly to the comparison of MadDM results on direct detection to the results in 
literature. 

The diagrams contributing to DM-nucleon cross section are displayed in Fig. The Higgs- 
exchange diagram contributes to SI DM-nucleon cross section while the other two diagrams in¬ 
volving KK quarks contribute both for SI and SD DM-nucleon cross sections. We compared the 
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Figure 2. Spin-independent elastic scattering cross section of scalar DM in the simplified SM model scenario 
with proton and neutron for <5 = 0.1 (top). Spin-dependent elastic scattering cross section of Dirac fermion 
DM with proton and neutron for = 0.1 (bottom). 



Figure 3. Feynman diagrams for elastic scattering of KK photon with a quark. 


MadDM results against those using the private code used in Refs. |46l I49j . The results, shown in 
Fig. show perfect agreement between MadDM and the private code over a wide range of dark 
matter masses and for several values of A. The DM-neutron cross sections are also in excellent 
agreement. 
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Figure 4. Spin-independent (left panel) and spin-dependent (right panel) DM-proton cross section for 
MUED. The blue curves correspond to theoretical values coming from Kaluza-Klein dark matter literature 
[1S1I35]. The green points are the MadDM data. We show the results for three values of the A parameter. 


3. Minimal Supersymmetric Standard Model 


Finally, we also validated MadDM for SPSla benchmark point in Minimal Supersymmetric 
Standard Model (MSSM) [50] . MSSM neutralino dark matter is similar to KK photon dark matter 
in that the processes relevant for direct detection typically proceed via squark exchange in the s- 
and u-channels and Higgs exchange in the t-channel (see Fig. [^. In Table IV, we compare MadDM 
with MicrOMEGAs for DM-nucleon cross-section for several SPS points. The ratio of the two tools 
cross-sections indicates that MadDM agrees with MicrOMEGAs to ~ 1% level, with the exception 
of the SPSS point. We have checked that the statement is also true if we vary the neutralino mass 
around the SPSla point (m^o = 96.68GeV). 


We note that it is challenging to compare other MSSM parameter points between MicrOMEGAs 
and MadDM. The default implementation of the MSSM model in MicrOMEGAs contains numerous 
model specihc optimizations and improvements which we are not included in MadDM. Hence, the 
only way we could compare the two codes in an “apples to apples” fashion in the context of MSSM 
was to use the MSSM implementation provided on the EeynRules website. Despite numerous 
technical difficulties, we were able to generate MSSM models for several specific SPS points, by 
first pre-processing the SLHA cards from SUSYHIT to ht the format required by FeynRules and 
then using those parameter sets to output both UFO and CalcHEP models. For the purpose of 
comparison, we manually turned off running of the QCD coupling in MicrOMEGAs by setting the 
variable qcdNL0=l in directDet. c and turned off the contributions from higher twist operators 
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Figure 5. Nuclear recoil energy (left panel) and angular (right panel) distributions for spin-independent 
interactions for different materials, assuming a 100 kg detector measuring events over one year for a DM 
mass of 100 GeV and DM-nucleon cross-section of 1 x 10“® pb. 

by setting Twist20n=0, as these features are not built into MadDM. 


B. Recoil Rates 


Upon validating the DM-nucleon scattering cross section results from MadDM, we proceed to 
the recoil rates for DM scattering off a target nucleus. We begin with a simple, model independent 
validation of the recoil rate calculation, where we simply assume that the DM-nucleon cross section 
a^n = 10® pb, chosen for the purpose of comparison with the results from Ref. [20]. To reproduce 
the SI recoil rates as a function of energy/angle as in Ref. [20|, we employ the differential recoil 
spectrum of Eq. (24), integrated over time and angle/energy. Fig. shows the spin-independent 
recoil rates as a function of recoil energy (left) and recoil angle (right). We find that both distribu¬ 
tions are in a very good agreement with the results found in Ref. [2^ , over a wide range of target 
materials. 

As a next validation, we check the recoil rates in UED model following the procedure described 
in Ref. [l8|, which shows sum of SI and SD recoil rates. The spin-dependent recoil rates are 
sensitive to numerical values of various quantities such as magnetic moments and parametrization 
of form factors. We use those values quoted in the references that are cited in Ref. |48j . Fig. 
shows nuclear recoil energy distributions as a function of recoil energy for Xenon, Germanium and 
Nal. KK photon mass is chosen to be 1000 GeV with the DM-nucleon scattering cross-sections for 
both spin-dependent and spin-independent for A = 15% as illustrated in Fig.|^ Despite the minor 
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cr type 

MadDM [pb] 

micrOMEGAs [pb] 

% difference 


SPSla 


SI (proton) 

2.16E-10 

2.I8E-10 

0.9 

SI (neutron) 

2.15E-10 

2.I7E-10 

0.5 

SD (proton) 

6.53E-6 

6.56E-6 

0.5 

SD (neutron) 

8.79E-6 

8.79E-6 

0.0 


SPSlb 


SI (proton) 

3.87E-1I 

3.89E-11 

0.5 

SI (neutron) 

3.85E-1I 

3.87E-11 

0.4 

SD (proton) 

I.20E-6 

I.2IE-6 

0.5 

SD (neutron) 

I.47E-6 

I.47E-6 

0.0 


SPS2 


SI (proton) 

7.35E-1I 

7.39E-11 

0.5 

SI (neutron) 

7.35E-1I 

7.39E-11 

0.5 

SD (proton) 

9.86E-6 

9.89E-6 

0.3 

SD (neutron) 

8.05E-6 

8.05E-6 

0.0 


SPSS 


SI (proton) 

5.55E-1I 

6.02E-11 

0.4 

SI (neutron) 

5.53E-1I 

5.98E-1I 

0.1 

SD (proton) 

I.63E-6 

I.64E-6 

0.6 

SD (neutron) 

2.08E-6 

2.08E-6 

0.0 


SPS4 


SI (proton) 

2.04E-10 

2.07E-10 

1.3 

SI (neutron) 

2.03E-10 

2.05E-10 

1.0 

SD (proton) 

7.56E-6 

7.58E-6 

0.3 

SD (neutron) 

7.82E-6 

7.82E-6 

0.0 


SPSS 


SI (proton) 

6.97E-11 

4.33E-11 

4.6 

SI (neutron) 

6.95E-11 

4.31E-11 

4.0 

SD (proton) 

4.05E-8 

4.09E-8 

1.0 

SD (neutron) 

4.77E-7 

4.77E-7 

0.0 


SPS9 


SI (proton) 

5.01E-11 

5.03E-11 

0.3 

SI (neutron) 

5.00E-11 

4.98E-11 

0.3 

SD (proton) 

l.llE-6 

1.12E-6 

0.6 

SD (neutron) 

1.72E-6 

1.72E-6 

0.0 


Table IV. MadDM and MicrOMEGAs comparison for DM-nucleon cross sections in MSSM. 
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Figure 6 . Nuclear recoil energy distributions as a function of recoil energy for Xenon, Germanium and Nal. 
These are the rates for the UED model assuming a mass of 1000 GeV for the DM ( 71 ) with the DM-nucleon 
scattering cross-sections for both spin-dependent and spin-independent for A = 15% as illustrated in Fig. 



MicrOMEGAs 

MadDM 


5.41 X 10-5 

5.42 X 10-5 

73Ge 

2.99 X 10-5 

3.07 X 10-5 

127j 

5.28 X 10-5 

5.30 X 10-5 

23Na 

3.09 X 10-5 

3.12 X 10-5 


Table V. Total number of events for a Higgs portal dark matter model with 100 GeV mass. We compare 
the number of events as obtained from MicrOMEGAS and MadDM for a 1kg detector measuring events 
over a day. 

differences in the spin-dependent recoil rates, we find that Fig. shows a good agreement between 
MadDM and the calculations in Ref. |48| . 

Furthermore, we compare the differential energy rates obtained from MicrOMEGAs with those 
obtained from MadDM for a Higgs portal scalar dark matter model for 4 different materials: Xenon, 
Iodine, Germanium and Sodium. We consider only the most abundant isotopes of these materials 
as is assumed in MicrOMEGAs. SD cross section vanishes in this case since DM is a scalar. Eig. 

illustrates the comparison between MicrOMEGAs and MadDM for the differential energy rates 
as a function of recoil energy in keVnr, while Table [V| shows several numerical comparisons of the 
total expected number of events for a 1 kg-day normalization. We find an excellent agreement 
between the two codes over a wide range of recoil energies and target materials. 
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Figure 7. Comparison of recoil energy distribution for I and Ge (left panel) and for Xe and Na (right panel) 
for a Higgs portal scalar dark matter model. 

C. LUX exclusion bound 

As a final validation of the MadDM code, we attempt to reproduce the exclusion limit on the 
DM-nucleon cross-section as a function of DM mass similar to the LUX 2013 experimental results 
|18| . For this purpose, we assume the efficiency function of nuclear recoils displayed in the black 
curve of Fig. 1 in Ref. ^5] . 

Fig. [^shows the results from MadDM for different energy threshold cuts as compared to the data 
reported by the LUX experiment. We present the contours assuming 2.3 events, coinciding with 
the number of events at 90% conhdence as required by the Feldman-Cousins confidence intervals. 
We find a good match between the LUX data and limits from MadDM. As we could not obtain 
information on what value of the energy threshold cut is used in the LUX limit, we considered 
different values of threshold cuts. As illustrated in Fig. the threshold cuts only impact the lower 
mass side since a higher threshold cut reduces the statistics for a lower DM mass. 

Note that, as described in the previous sections, the exclusion curves in Fig. [^can be obtained 
in MadDM by using LUX_Exclusion routine found in the test routines part in maddm.f. The 
routine multiplies by the efficiency obtained from Ref. [23], which is is then weighted by a 
50 % acceptance rate for nuclear recoils as stated in the LUX analysis. From the recoil spectrum 
weighted by the efficiency and the acceptance rate, the function then calculates the total number 
of expected events. The default value for the detector efficiency is 100%, and can be easily replaced 
by a user defined function. 
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Figure 8. 90% confidence limits on the spin-independent DM-nucleon scattering cross-section (in picobarns) 
for an unsmeared energy distribution (left panel) and the smeared distribution with A = 1 (right panel). 
Limits are obtained from MadDM for 2 keV (black solid) and 4 keV (red dashed) and LUX limits are shown 
in blue curve with circular data points. 


V. CONCLUSIONS 

The identity of dark matter is one of the most profound mysteries in particle physics, astro¬ 
physics, and cosmology. Recent data from gamma rays, supernovae luminosities, cosmic microwave 
anisotropies, and galactic rotation curves all point consistently to the existence of dark matter with 
~ 5 times more abundance compared to ordinary matter. At the same time, all known particles are 
excluded as possible dark matter candidates, making the dark matter problem perhaps the most 
pressing motivation for physics beyond the Standard Model. Little is currently known about the 
mass scale of dark matter, suggesting that discovery and characterization of DM will likely require 
a synergistic approach including well-balanced programs in direct detection, indirect detection, 
particle colliders and astrophysical probes. 

In order to efficiently combine results from various dark matter searches sparked a demand 
for a new generation of numerical tools. MadDM is an on-going effort to bridge DM collider 
phenomenology with astro-physics and cosmology of DM, with the ultimate goal to provide an “all 
in one” dark matter phenomenology package which can be easily incorporated into the future dark 
matter searches at the LHC. 

In our current work, we presented MadDM 2.0, which includes direct detection of dark matter 
in a generic UFO model. The code computes the total DM-nucleus scattering rate, recoil energy 
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spectrum and angular distributions of dark matter elastic scattering off nuclei for various target 
materials, including user-provided detector resolution and detector efficiency. MadDM v.2.0 follows 
the “easy to use, easy to install” philosophy of the previous version, featnring great flexibility so that 
nsers can replace the existing modnles with their own implementation or link MadDM functionality 
with other existing tools. We have performed detailed valnation of MadDM against MicrOMEGAs 
and private codes for various benchmark dark matter models and found excellent agreements with 
existing tools and literature. 

In addition to relic abundance and direct detection, the fnture renditions of the MadDM code 
will also provide ability to calculate cosmic-ray fluxes from DM annihilation, including gamma 
rays, cosmic electrons/positrons and protons/anti-protons. Fnrthermore, recent improvements in 
aMC@NLO framework allow for compntations of amplitndes for loop induced diagrams , giving 
ns an opportunity to create the first automated DM tool which can calculate relic density, direct 
detection and indirect-detection signals originating in loop-indnced processes. 
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Appendix A: Structure of MadDM v.2.0 

MadDM v.2.0 follows the directory and hie structure of MadDM v.1.0. The Python module is 
arranged into the following simple structure: 

• EffOperators: A folder containing UFO hies for effective operators, which are used to 
extract the low energy coefficients of the XQ. XQ matrix elements. The UFO hies were 


generated using the FeynRnles 2.0 convention (see Section IIIB for more details.). 


ExpData: all experimental data used by MadDM is stored in this folder. Gurrently, only the 
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LUX exclusion bound for the spin-independent DM-nucleon cross section is included. Future 
versions of MadDM will include a more elaborate database of experimental bounds. 

• All user defined projects are automatically output into the Projects folder. All numerical 
calculations and result output is written into the Projects folder. 

• Templates contains necessary FORTRAN routines and template files for relic density and direct 
detection. The structure of the Template folder serves as a skeleton for a user defined Project. 

• MGoutput. py contains a simplified version of the MadGraph output routines, which write 
the matrix elements into FORTRAN files. The routines within MGoutput.py will write the 
matrix.! files, into the Projects folder, but will omit writing the additional files which are 
not necessary for MadDM to run, but are ordinarily written by MadGraph output routines. 

• darkmatter .py contains the essential part of MadDM Python module. The routines within 
this library contain class definitions for the darkmatter objects, as well as numerous func¬ 
tions which find dark matter candidates and generate the MadDM FORTRAN module. 

• A self contained example of a MadDM calculation without the use of the user interface can 
be found in example.py. The example calculation goes through the process of importing 
a user defined UFO model, finding the DM candidate, generating the FORTRAN module and 
calculating the DM-nucleon cross section and relic density. 

• init . py contains the definitions and functions used by the MadDM user interface. 

• maddm. py is the main Python executable. It will initiate the user interface defined in init. py 
and lead the user through the steps of DM relic density and calculations relevant for DM- 
nucleon scattering. 

• param_scan_def ault .py is a pre-defined parameter scanning Python script. The file con¬ 
tains a skeleton of the Python code which will scan over DM model parameters and calculate 
relevant physical quantities. In order to run the parameter scanning script, a minimum of 
user intervention is required. All places where user input is necessary are marked, and lim¬ 
ited to the choice of DM model parameters, number of parameter to scan over and definitions 
of parameter ranges. 

The FORTRAN module in each Projects folder contains the following 7 subdirectories. 
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• Cards: A directory which contains param_card.dat. The structure and placement of the 
parameter card follows the standard formats of the MadGraph5_aMC@NLO framework. 

• Source directory contains FORTRAN code for the DM model, as well as necessary HELAS 
libraries. 

• include directory contains various include files for the FORTRAN code. The maddm_card. inc 
which contains user defined parameters of the DM relic density and direct detection calcu¬ 
lations is located in the include folder. Note that any changes in the maddin_card. inc file 
will only be reflected upon the re-compilation of the MadDM FORTRAN module. 

• Amplitudes for relic abundance and direct detection are saved under matrix_elenients. 

• All results (relic density, cross sections and distributions) are found in output directory. 

• src directory contains the MadDM FORTRAN code. The main program is located in the 
maddm.f file. 

All the other files in the main directory are executable and perform the following tasks. 

• maddm. x is the main FORTRAN executable. Note that all output of the maddm. x will be written 
into the output/maddm. out file by default. 

• make_plots is a Python script containing calls to plotting MadDM plotting routines. All 
calls are commented out by default. Note that it is necessary to install both Numpy and 
Matplotlib in order for plotting functionality to work in MadDM. 

• plotting.py is a library of user friendly plotting functions. See make_plots for examples 
of how to use plotting functionality of MadDM. Function definitions in plotting.py use 
libraries from both Numpy and Matplotlib. 

Appendix B: Description of MadDM v.2.0 rontines 
1. Python Module 

We upgraded the Python module of MadDM to allow users to select which calculation they 


would like to perform. For this purpose we introduced three switches, which are members of 
darkmatter class: 
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• <darkmatter> _do_relic_density: True if the user wants to calculate dark matter relic 
density, False otherwise. 

• <darkmatter> _do_direct_detectioii: True if the user wants to calculate spin-dependent 
and spin-independent cross sections for dark matter scattering off protons and neutrons, 
False otherwise. 

• <darkmatter> _do_directional_detection: True if the user wants to calculate recoil rates 
(as a function of energy and angle) assuming detectors made of various materials. Please 
note that this flag can not be set to True if _do_direct_detection = False. 

In order to reflect the changes in the code, we have changed the name of the function which 
executes the numerical module of MadDM from CalculateRelicDensity() to CalculateO. The 
function now takes the output file name, output_file as a parameter, specifying the path of 
the file where all the calculation results will be written. The default value for the output file 
is output/maddm.out. Changing the output file path is particularly useful if the user wishes to 
parallelize a parameter scan, as it avoids the problem of many jobs attempting to read/write the 
same file. 

The only two new relevant Python functions in the darkmatter class are: 

• darkmatter: [excluded, excluded_omegah2, excluded_dd]: 

is_excluded( omega_min, omega_max , si_limit_x, si_limit_y ): This function 

compares the results of a DM calculation obtained with the CalculateO function with 
the existing bounds on relic density and spin-independent DM-nucleon scattering cross sec¬ 
tions. The parameters omega_min and omega_max represent the minimum and maximum 
for which the parameter point is not ruled out while the si_limit_x and si_limit_y are 
one dimensional arrays of numbers representing the x and y values of the spin-independent 
cross section exclusion limit. The default limit implemented in MadDM is the LUX bound, 
which can be found in a text file inside the ExpData subdirectory of MadDM. The function 
assumes that the DM-nucleon cross section values are given in pb. The function 
returns an array of three boolean variables, excluded signals whether the model point 
is excluded overall, while excluded_omegah2 and excluded_dd signal whether the model 
point is excluded by the individual constraint. This function should be called only after 
CalculateO has been executed. 
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• darkmatter: None: GenerateDiagramsDirDetect(): The routine generates all XQ. XQ 
diagrams, where q = u,d, s, c, b, t quarks, as well as XQ XQ diagrams and extracts the 
interference terms relevant for computing the low energy effective operators for x-nucleon 
scattering. The function takes no values and returns no valnes. 

Simply calling GenerateDiagramsDirDetect () will automatically set _do_direct_detection 
to True and result in a calculation which will go as far as y-nucleon spin-independent and spin- 
dependent cross sections. The user can proceed to calculate the differential recoil rates and simulate 
dark matter events scattering off nuclei by also setting the flag _do_directional_detection = True 


2. FORTRAN Module 


In the FORTRAN module, we added a library of function definitions in direct_detection. f and 
directional_detection.f files inside the src directory, which calculate the direct detection nu¬ 
cleon cross sections, differential recoil rates and other relevant qnantities. The relevant function 
definitions are: 


• sigma_nucleon(nucleon) : The function calcnlates the numerical value of a cross section for 
dark matter particle scattering off a nucleon at 0 momentum transfer. nucleon=l,0 refers 
to a proton/neutron respectively. The function reads in the nucleon form factors as defined 
in maddm_card. inc (see below for more detail). 


FHelm(Er, A) : This function calculates the Helm form factor as given in Eq. (25). The 
function reads the recoil energy in keVnr (Er), and the mass nnmber of the target atom (A) 
as input parameters. 


• FWS(Er, A) : This function calculates the Wood-Saxon form factor. The function reads the 
recoil energy in keVnr (Er), and the mass number of the target atom (A) as input parameters. 

• FWSmass(m_DM, A): This function calculates the Wood-Saxon form factor. The function 
reads the WIMP mass in GeV (m_DM), and the mass number of the target atom (A) as input 
parameters. 


coefficients (mater, coeff_s00, coeff_s01, coeff_sll) : This routine outputs the coef¬ 
ficients for the polynomial fits to the nuclear structure functions Sij for the spin-dependent 


form factor. These are referenced in table III The input is the material mater which is 


chosen in maddm_card. inc. 
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• Spiii_matrix (mater, J, avg_sp, avg_sn) : Provides the spin J and nuclear magnetic mo¬ 
ments avg_sp and avg_sn for the different target nuclei, mater is the material chosen in 
maddm_card.inc. 

• Structure_Func(Er, A, mater, sOO, sOl, sll) : Uses the coefficients routine to cal¬ 
culate the structure functions (sOO, sOl and sll) that are used in the spin-dependent nuclear 
form factor. The input is the recoil energy Er in keVnr, the atomic mass A and the target 
material mater. 


FormSD(Er, A, mater) : Calculates the spin-dependent form factor as a function of recoil 


energy in keVnr from the structure functions in Structure_Func according to Eq. (26). 


FormSDmass(m_DM, A, mater) : Calculates the spin-dependent form factor as a function of 


WIMP mass in GeV from the structure functions in Structure_Func according to Eq. (26). 


• VE(days): This function calculates the velocity of the earth (in km/s) relative to the galactic 
center as a function of time in days according to appendix B in mi- The function takes 
into account the proper motion of the solar system with respect to the surrounding stars, 
the motion of the Earth around the Sun and provides galactic coordinates of the Earth as 
it moves each day on its orbit. The input parameter days represents the day of the year 
relative to December 31*^* of the calendar year. 


d2RdEdcostheta 

(Er, costheta, sigmawnOSI, sigmawpOSI, sigmawnOSD, sigmawpOSD, M_dm, ve) : 
This routine calculates the recoil rates as a function of energy and angle on a particular 


day of measurement, based on Eq. (24). The maddm_card, inc file contains numerical val¬ 
ues of necessary quantities such as the target material, size of the detector (in kg), the 
velocity of the earth (calculated inside VE) , the most probable DM velocity uq, the es¬ 
cape velocity of the DM from the galactic Halo Vesc, the minimum DM velocity required 
for recoil Vmin all in units of km/s, and the local halo density poj in GeVc“^cm“^. The 
variables sigmawpOSI and sigmawnOSI represent the spin-independent cross sections, while 
sigmawpOSD and sigmawnOSD are the spin-dependent cross-sections for DM scattering off 
protons and neutrons respectively in units of pb. M_dm and ve are the dark matter mass in 
GeV and the Earth’s velocity relative to the Galactic center respectively. 
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• target_material(material, A, Ab, Z, n) : This subroutine provides the double differ¬ 
ential module d2RdEdcostheta with the correct target material the user has chosen inside 
maddm_card. inc. The input is an integer material provided through a common block from 
the maddm_card. inc file, while the output is an array of mass nnmbers (A) for the target 
atoms, their respective abundances (Ab), the atomic number (Z) as well as the length of the 
mass number arrays (n). These are listed in Table [TT} For instance, consider that the target 
material has been set to 1 (Xenon, which is default) in the maddm_card.inc file. Then 
target_material will output an array of 9 elements for the mass number of the 9 most 
stable Xenon isotopes and their corresponding abundance information. It will also provide 
Z = 54. 

• GaussmearAngle(uaevents, smaevents) : This function smears a histogram of theoretical 
distributions (uaevents) nsing a Gaussian distribution. The angular resolution is set inside 
maddm_card. inc by changing the sig_theta parameter m- The default angular resolution 
is 30°. 

• GaussmearEnergy (uevents, smevents) : This module applies energy smearing to the dis- 
tribntion of angle smeared events (uevents) using a Gaussian resolution function with un¬ 
certainty aE = X\/Ee,, where Er is the recoil energy in keV and A is a constant which 
determines the level of smearing. The default value is A = 1 inside maddm_card. inc. The 
smearing function can convert low energy recoil events into negative energies. An energy 
threshold cut alleviates this problem. 

• Theory_Simulation 

(mass, sigma_wnSI, sigma_wpSI, sigma_wnSD, sigma_wpSD, Theory): This subrou¬ 
tine uses the dark matter mass in GeV (mass) and DM-Nucleon cross-sections for the 
spin-independent ( sigma_wpSI for proton, sigma_wnSI for neutron ) and spin-dependent 
( sigma_wpSD for proton, sigma_wnSD for nentron ) contributions to simulate the theo¬ 
retical double differential distribution (with respect to recoil energy and recoil angle) and 
outputs the result as an array named Theory representing the distribution of unsmeared 
recoil events. In this routine, the double differential rate is approximated at the center of 
each energy-angle-day bin. 

• ef f (x) : Function that interpolates the data points from LUX for the nuclear recoil efficiency. 
The input x is the recoil energy (in keV) and the output is the detector efficiency at that 
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energy. 

• efficiency(x, flag) : Provides the DM detection efficiency as a function of recoil energy. 
The variable x represents the recoil energy in keV and flag is an integer variable which, when 
chosen as 0, gives a constant 100% detection efficiency. If flag is 1, the efficiency is 50% 
and if flag is 2, the efficiency is that obtained from the function eff, which is the nuclear 
recoil detection efficiency from the LUX experiment in Fig. 1 of Ref. [23]. For a conservative 
estimate we extend and interpolate the efficiency to 100 keV. New detector efficiency may 
be implemented by user in place of this function definition. 

• directional_detection 

(mass, sigma_wnSI, sigma_wpSI, sigma_wnSD, sigma_wpSD, flag, N_events): This 
routine calculates the differential distributions of dark matter recoil rates given the spin- 
independent ( sigma_wpSI for proton, sigma_wnSI for neutron) and spin-dependent 
( sigma_wpSD for proton, sigma_wnSD for neutron) cross sections in pb as well as the 
DM mass in GeV (mass). The subroutine calls the theory_simulation routine which 
performs the event distribution simulation. If the user specified smearing = . true. in 
the maddm_card. inc file the function will incorporate detector smearing into the resulting 
distributions. The function calculates the double differential distribution (as a function 
of energy and angle), and integrates over angle and time to obtain the differential energy 
distribution (^). The energy distribution is weighted with the detector efficiency function 
efficiency. Furthermore, the function also integrates over energy and day to obtain the 
differential angular distribution (y^^ ) as well as the total rate {^ = R) after integrating 
over energy and angle. All the distributions are written to separate files in the user defined 
project output directory. The total number of events for a 1 ton detector over the course of 
1 year is stored in the variable N_events. The integer variable flag is used to choose the 
type of efficiency to be used. 

In addition to functions which calculate quantities relevant for DM-nucleon scattering, we added 
a number of important variables to the FORTRAN module, in the maddm_card. inc file inside the 
include directory of the project folder]^: 

• do_relic_density: A logical flag which determines whether a relic density calculation 
should be performed. It can be either .true, or .false. If users have not gener- 

® Please note that one must recompile the numerical code by typing make in the project folder in order for any 
changes to the maddm_card.inc hie to take effect properly. 
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ated the FORTRAN module with the relic density calculation [i.e. there is no call to 
GenerateDiagrams0 in the Python module), setting this flag to .true, will result in 
an error. 

• do_direct_detection: A logical flag which determines whether a direct detection calcu¬ 
lation should be performed. It can be either .true, or .false. If users have not gen¬ 
erated the FORTRAN module with the direct detection calculation (i.e. there is no call to 
GenerateDiagramsDirDetect0 in the Python module), setting this flag to .true, will 
result in an error. 

• do_directional_detection: A logical flag which determines whether to calculate differ¬ 
ential distributions for recoil rates (both energy and angle dependence). It can be either 
.true, or .false. If users have not generated the FORTRAN module with the direct de¬ 
tection calculation (i.e. there is no call to GenerateDiagramsDirDetect() in the Python 
module), setting this flag to .true, will result in an error. In addition, one can not set to 
do_directional_detection = .true, ifdo_direct_detection = .false.. 

• SPu. . . SPg, SNu. . . SNg: Values for the scalar nucleon form factors for the proton (SPx) 
and neutron (SNx). The nucleon constituent {u,d,...g) is labeled by the last letter in the 
variable name. 

• VPu, VPd, VNu, VNd: Values for the vector nucleon form factors for the proton (VPx) and 
neutron (VNx). The nucleon constituent {u, d, ...g) is labeled by the last letter in the variable 
name. 

• AVPu. . . AVPs, AVNu. . . AVNs: Values for the axial-vector nucleon form factors for the 
proton (AVPx) and neutron (AVNx). The nucleon constituent {u,d, ...g) is labeled by the last 
letter in the variable name. 

• SigPu. . . SigPs, SigNu. . . SigNs: Values for the tensor nucleon form factors for the 
proton (SigPx) and neutron (SigNx). The nucleon constituent {u,d,...g) is labeled by the 
last letter in the variable name. 

• material: Variable which selects the detector material. It can take values 1-13 for the 
following implemented materials; 1 - Xe, 2 - Ge, 3 - Si, 4 - Ar, 5 - Ne, 6 - Na, 7 - I, 8 - C, 
9 - F, 10 - S. For compound materials set mater_comp as .true, in maddm_card.inc, 11 - 
Nal, 12 - CF 4 and 13 - CS 2 . 
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• vMP: Most probable speed for dark matter particles in the galactic halo in km/s, assuming 
the Maxwell-Boltzmann velocity distribution. The default value of vMP is vq = 220 km/s. 

• vescape: Escape velocity for dark matter in the galactic halo in km/s, the default value is 
Uesc = 650 km/s. 

• rhoDM: Mass density of dark matter in the halo in GeVc“^cm“^, default value is 0.3 

GeVc“^cm“^. 

• detector_size: Mass of the dark matter detector target in kg, the default value is 1000 kg. 

• En_threshold: The recoil energy threshold cut (in keVnr). 

• lambd : Energy resolution parameter defined in aE = Av/E with A = 1 as default. 

• sig_theta: Angular resolution of the detector (in degrees) with default value of ag = 30°. 

• smearing: A logical flag which determines whether the results should be smeared by de¬ 
tector effects. If set to . true. it will smear both the recoil energy and angle. 

Eor the purpose of dark matter scattering event simulations, we also incorporated a number of 
parameters which are used as inputs for the nuclear recoil scattering simulation: 

• En_min, En_max: Minimum and maximum recoil energy in keVnr. 

• cos_min, cos_max: Minimum and maximum cosine of the recoil angle. 

• day_min, day_max: Starting and ending day of the experiment relative to December 31®*. 

• Energy_bins, cos_theta_bins, day_bins: The number of bins in histograms for scat¬ 
tering rate as a function of recoil energy, cosine angle and day of experiment respectively. 

3. Test routines 

Eor debugging purposes we also provide a set of test routines and functions which can be used 
in maddm.f and are commented out by default. Simply un-commenting the calls and re-compiling 
the FORTRAN module will make the test functions accessible. The list of test function includes the 
following: 
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d2RdEdcos_test 

(Er, costheta, sigmawnOSI, sigmawpOSI, sigmawnOSD, sigmawpOSD, M_dm, day): 
This subroutine uses the d2RdEdcostheta routine described in the previous section to print 
out all the ingredients that go into the double differential function, for consistency checking. 
The input variables are defined in d2RdEdcostheta above. The output is printed on the 
screen and includes: DM mass, the spin-dependent and spin-independent cross-sections, the 
target material used as well as its mass numbers and the isotope abundances, the detector 


size and other variables that factor into Eq. (24). 


• V_E_test (day, nuni_array) : This routine provides the velocity of the Earth through the 
galactic halo on a specific day relative to December 31®*. nuni_array is an integer flag. When 
set to 0, V_E_test prints one number representing the velocity on the day specified by day. 
However when num_array is set to 1, the routine outputs an array of days and velocities (in 
the galactic frame) with the minimum and the maximum days set in the niaddm_card. inc 
file. The output file is written in output/VE_test. dat. 


• Form_test(Er, num_array) : Outputs the spin-dependent and spin-independent form fac¬ 
tors for the most abundant elements. It uses the subroutine target_MA, which is exactly 
the same as the routine target_niaterial above, but only provides the most abundant 
elements. When num_array is 0, the output is the form factor at a specific input recoil 
energy. When num_array is 1, the output is an array of the form factor squared in the file 
output/f ormf ac_test. dat. The first column is the recoil energy, the second is the spin- 
independent form factor squared while the third is the spin-dependent form factor squared. 


• Form_testmass (ni_DM, num_array): Outputs the spin-dependent and spin-independent 
form factors for the most abundant elements. It uses the subroutine target_]yiA. When 
nuni_array is 0, the output is the form factor at a specific input WIMP mass in GeV. 
When num_array is 1, the output is an array of the form factor squared in the file 
output/f ormf acniass_test. dat. The first column is the WIMP mass, the second is the 
spin-independent form factor squared while the third is the spin-dependent form factor 
squared. 


LUX_Exclusion(xmim, xmax, ymin, ymax, step) : Scans the DM-nucleon vs Mass param¬ 
eter space to provide the LUX exclusion limit. This is done by using the flag = 2 in the 
directional_detection routine, which ensures that we use the efficiency function obtained 
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from the LUX 2013 results [23]. xmin represents the minimum mass, xmax, the maximum 
mass, ymin is the minimum cross-section and ymax is the maximum cross-section in the 
scanning parameter range, nstep is the number of scanning points for each direction. The 
scanning is done in log space, the linear scanning is commented out, for linear scanning, one 
has just to uncomment the Lux_exclusion line in the test routine part of maddm.f 
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